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CTN , We study by means of numerical simulations the model of dissipative particle dynamics with 

Q j ' energy conservation for the simple case of thermal conduction. It is shown that the model displays 

correct equilibrium fluctuations and reproduces Fourier law. The connection between "mesoscopic 

coarse-graining" and "resolution" is clarified. 
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I. INTRODUCTION 
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Mesoscopic approaches to the study of complex fluids are receiving a great deal of attention in view of their potential 
r^ , to address problems with disparate time-scales. The method of dissipative particle dynamics introduced originally as 

' an hybrid between lattice gas and molecular dynamics simulations |l| , has proven to be a flexible and powerful tool in 
^ [ the study of complex fluids as colloidal suspensions f^-^ , porous flow [y , polymer suspensions |^ , and multicomponent 
S ' flows §. 

1 K I DPD models a fluid as a collection of point particles that are interpreted as representing the center of mass 
C^ . of a mesoscopic cluster of the atoms or molecules that constitute the fluid. By formulating physically motivated 
c/3 interactions (that is, dissipative forces that conserve momentum [L7]) between these "lumps" of fluid, it is assured 

that the macroscopic behavior of the system is hydrodynamic |B,3|. Most importantly, a random interaction is 
introduced in order to account for thermal fluctuations that occur at the mesoscopic level and are responsible for 
Brownian effects appearing in complex fluids. 

One of the drawbacks of classic DPD is that the total energy of the system is not conserved because the forces 



^ ' are velocity dependent. This has been recently remedied simultaneously and independently in Refs. [^ 11 1, where an 

Q internal energy variable and a temperature is introduced for each particle. The mechanical energy that is dissipated 

^ due to the velocity dependent forces is transformed into internal energy of the particles. This viscous heating process is 

'~~' , accompanied by a thermal conduction process where exchange of internal energy between neighboring particles occurs 

, ' when differences of temperatures exist. The new DPD algorithm with energy conservation opens up the possibility of 

^ , studying not only rheological aspects of complex fluids but also thermal issues. It is therefore important to validate 

^-H ' the technique in simple situations for which the dynamical aspects are well-known, either analytically or by means of 

CNJ other numerical simulations techniques. 

^~^ As a first step towards the understanding of the behavior of DPD with energy conservation we focus in this paper on 

^N ■ the simplest problem of conduction in a quiescent fluid (or a solid) . This is a particular case of the model introduced in 

ir . Refs. 10,1^]. The particles are at rest and located randomly and represent portions of material that can interchange 

^■^ ' energy. This simple model can be regarded as a way of solving a fluctuating heat conduction equation. In section 

■*-«. i 2 we review the model that was introduced in Refs. |lfl,0], for the case that only heat conduction takes place. In 

j^ ' section 3 we further specify the model functions. In section 4 we present simulation results that show the validity 
of the simulation method by comparing the numerical results with the theoretical results. Finally, we end up with a 
section of conclusions. 



II. THE MODEL 



. !^ ' We have introduced in Ref. |T^] a mesoscopic model for describing a fluctuating viscous and thermally conducting 

^^ fluid. Here we apply this model to a quiescent fluid (or solid) within a box. The system is represented by a set of 

^H N particles of fixed positions r^ which are distributed at random within the box. These particles are understood 

. 5t 1 as thermodynamic subsystems of the whole system, that is, small portions of material that have a sufficiently large 

number of degrees of freedom. To each particle we associate an internal energy variable Ci an entropy function s{ei), 

and a temperature T^ = dsi/dci. Each particle represents a mesoscopic portion of the material and its internal 

energy content is subject to thermal fluctuations, which in a continuum theory would be described by a random heat 
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flux. [Q Differences in temperatures between neighboring particles wifl cause a transfer of internal energy. Therefore, 
the following stochastic differential equation for the internal energy of each particle is postulated 

j \ * J / j 

The first term in the r.h.s. is deterministic and specifies that a temperature difference causes exchange of energy. The 
second term is stochastic and takes into account thermally induced fluctuations in each particle caused by a random 
interchange of energy between particles. The notation is as follows, ks is Boltzmann's constant, r^ — jr^ — r^l is the 
distance between particles i,j and the dimensionless weight function tu{r) determines the range of influence between 
particles. The parameter Kij governs the overall amplitude of the deterministic and random parts. It depends in 
general on the state variables of particles i, j and is symmetric under particle interchange. The random "heat fiux" 
from particle i to j is expressed in terms of the increments of the Wiener process dWl^ . The increments of the Wiener 
process are antisymmetric under particle interchange dW^j = —dWj^ in such a way that the total internal energy of 
the system governed by Eqn. (^ is exactly conserved ^ ^^ e^ = 0. Because Kij might depend in general on ei,ej we 
need to provide a stochastic interpretation of Eqn. (|l]), as the noise is multiplicative. We choose Ito interpretation. 
It has been shown that, if dijKij = OR, then Eqn. (nl) is mathematically equivalent to a Fokker-Planck equation which 
has as equilibrium solution 
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where T'i^f.ek) is a function of the total internal energy which is determined by the initial distribution of total 
energy. Particular examples are given by the canonical pil and microcanonical ensembles 
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Here, the factors ft{Eo,N) and Z{f3,N) are obtained by normalization. 

III. PERFECT SOLID AND FOURIER LAW 

In this section wc further specify the model by providing the undetermined functions T(e), K(e + e') and a;(r) 
appearing in Eqn. (^. 

The simplest possible form for the equation of state T(e) is that of a perfect solid, which is a very good approximation 
for metals. The perfect solid has the following equation of state 

r(6) - -^ (4) 

where the heat capacity of the mesoscopic particles C^ is a constant independent of the energy. It is an extensive 
property that depends on the "size" of the particles. For mesoscopic particles, the dimensionless heat capacity 
a = Cy/ks ^ 1. The entropy s(e) is obtained by integrating the temperature and, apart from irrelevant constant 
additive terms, it is given by s(e) — Cy Ine. 

Now, let us assume the following functional form for k.^. 



^'Rtt^ =T^m+r.r (5) 



'^This implies that if t^ij depends only on ei,ej then, necessarily, it will be a function of e^ + e^. This can be proved by 
considering a Taylor expansion of a function of two variables and using the fact that, at each point, both partial derivatives 
coincide. 



where A is the average distance between particles (this is A = n~^l'^ and n is the average density number) and k has 
dimensions of a diffusion coefficient. In the last equality, Eqn. (Q) has been used. 

Substitution of Eqn. (|^) into the dynamical equation for the energy leads to the following equation for the temper- 
ature of a particle 
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The factor IxT- should be very close to 1 if the temperature differences between particles within an action sphere 
are small. We see, then, that for the perfect solid by a proper selection of Kij the form of the postulated SDE for the 
energy leads to a discrete version of the Fourier law of heat conduction. Note that the factor A^ in the first term of 
the left hand side of Eqn. (g) can be understood as the "lattice spacing" squared that would arise in a finite difference 
discretization of the Laplacian in the heat conduction equation. In a sense, Eqn. (|6|) can be understood as a finite 
difference discretization of the macroscopic heat equation on a random lattice which, in addition, has the correct form 
for the thermal fluctuations. 

We will assume that the weight function ^(r) has the form of the Lucy weight function |13[ , which is a smooth 
bell-shaped function with continuous derivatives. It is given by 

105 /, „ r \ / ^ r ^ 



where Tc is range of interaction between particle and s = re/ X is the overlapping coefficient which gives, essentially, 
the number of neighbors of a given particle. We have normalized the weight function as ||I| 

d'^rcoir) = - (8) 

n 



IV. SIMULATION RESULTS 

In this section we present some results obtained from numerical simulation of Eqn. (ph. The physical system 
is assumed to be in a three dimensional cubic box of edge length L. Initially the box is seeded with N mesoscopic 
particles located at random. The density is then n = N/Li^ and the typical interparticle distance is A = n~^^^. Periodic 
boundary conditions are assumed in the y, z directions and in the x direction we assume either periodic boundary 
conditions when considering equilibrium situations or we impose the temperature at a; = —L/2 and x = L/2. This is 
achieved by considering two extra layers of particles in each boundary in the x direction filled with particles that are 
at a constant temperature. These layers act, therefore, as thermal baths. They have a width as large as re in order 
that any particle within the system interacts with the same number of particles. 

The basic model parameters are the following ones: k, N, L, re, Cy, Tq. Here, Tq is a reference temperature, that 
of equilibrium or, in the case of nonequilibrium situations, that of the colder thermal bath. We work in reduced units 
such that L is the unit of length, L- /ii is unit of time, Tq is the unit of temperature, and ks is taken as the unit of 
entropy or heat capacity. 

We introduce the dimensionless heat capacity of mesoscopic particles a = Cy/ks and the dimensionless heat 
capacity of the total simulated sample, which is C^ — Na because heat capacity is an extensive property. On the 
other hand, Cl = L^Cy/ks, where c^ is the heat capacity per unit volume, which is a material constant. For example, 
Ci, = 3.44 X Kfj /Knn? for copper at room temperature. Therefore, by fixing Cl we are fixing the actual volume of 

the sample. For copper, a value Cl — 10® implies a submicron sample size of L = C^ x l.&\G~^'^ra ~ 0.1/i. In most 
of the simulations we use Cl = 10® and this fixes the value of a depending on the number of particles in the system. 
Note that by fixing Cl, large N implies small a, in accordance with the idea that the heat capacity of the mesoscopic 
particles is proportional to its typical size. We note that the overall intensity of the noise in (|g) is proportional to 
(/cb/C-u)^/^ = a~^/^ and, therefore, inversely proportional to the square root of the volume of the particles. Q The 
higher the number density of mesoscopic particles used to discretize a given sample, the smaller is the size of such 
particles and, therefore, the larger is the thermal noise. 

Eqns. (||) are solved with a conventional Eulcr method. For values of a < 10, the stochastic term produces from 
time to time a temperature which is negative for some particle. This has disastrous consequences because the factor 
accompanying (T^ — Tj) in the deterministic term of m) is negative and produces a flow of energy from the cold 



particles to the hot particles. When this happens the system becomes unstable. However, this never occurs for 
sufhciently large values of a as the ones we typically use and no instabilities are observed. 

In order to check that the code performs properly, we first run a simulation with fully periodic boundary conditions. 
Initially, all the particles are assumed to have the same temperature T — 1. After some equilibration period the system 



reaches equilibrium. The energy is exactly conserved (no round-off errors) because we use a Verlet list |15 , in such a 
way that what is gained by a particle is exactly what is lost by the rest of particles. 

In Fig. m we show the probability density distribution of the energy of a particle as obtained from simulation of 
Eqn. (yj) for two different values of a = C^/ks, the dimensionless heat capacity of the mesoscopic particles, for a 
system of A^ = 1000. The larger a the larger the average value of the energy and variance. However, the relative 
fluctuations are smaller, proportional to N~^/'^. Also shown in Fig. || are the canonical /can(e) and microcanonical 
/mic(e) theoretical results for the given values of a. The analytical expressions can be obtained from dsh and are given 
by 

/ca„(e) = ^i ^.. (/3e)"exp{-/3e} 
I (a + 1) 



where T{x) is the gamma function. The parameter f3 can be explicitly computed for a perfect solid by requiring 
that the average energy of a particle is Eq/N, where Eo is the total energy of the isolated system. The result is 
/3 = A'^(a + l)/Eo. In the thermodynamic limit the canonical and microcanonical distributions are identical. In 
practice, they are indistinguishable for values Na > 50. The perfect coincidence of the theoretical and simulation 
results in Fig. |l| gives confidence on the implementation of the numerical code used. 

Next, we compute the thermal diffusivity by a macroscopic measurement. Two heat baths are constructed with 
temperatures Tcom and Thot in such a way that a gradient (Thot — T^oidj/L is applied. A steady state is reached 
for sufficiently long times. Fig. shows the average of the temperature once the steady state has been reached for 
different number of particles for the case that Tcom — 1, Thot = 2 and the overlapping coefficient is s ~ 1.5. For 
A'' = 100 the spatial noise is considerable. Fluctuations in this graph are not due to statistical noise but to the spatial 
inhomogeneities that appear for small number of particles and they would smear out by averaging over different 
realizations of the initial spatial configuration of the particles. For larger number of particles a smooth linear profile 
is achieved, in accordance with our expectation that the discrete model reproduces Fourier law. The temperature 
field T{x) is obtained by binning the x axis and averaging the temperatures of the particles within each bin of volume 
AxL^. 

We have performed a series of simulations with different temperature gradients and have computed the macroscopic 
heat flux in the steady state. Two different ways of computing the macroscopic heat flux have been used. In the first 
case, the macroscopic heat flux is obtained by computing at each time step the energy gained by the cold bath (which 
is equal to the energy lost by the hot bath) and dividing it by the time step, this is AE/At. The second way follows 
from the microscopic definition of the heat flux as obtained from projection operator methods or from kinetic theory 

M 
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All transport fluxes are caused by coUisional transfer. As the particles do no move, there are no kinetic fluxes. This 
second form produces better results because it involves a very large number of pairs of particles and it is therefore 
adopted here. 

In Fig. we plot the macroscopic heat flux in terms of the applied temperature gradient. A linear dependence is 
obtained, which allows to extract the thermal diffusivity from the slope. 

When this experiment is performed for different number densities and different overlapping we can compare with 
the results from kinetic theory as developed in Ref . [|l^ . Fig. ^ shows the value of the macroscopic thermal diffusivity 
as a function of the overlapping s. Also shown is the kinetic theory prediction for the thermal diffusivity [M, reading, 

D = —k (11) 

We observe that the theoretical prediction improves for large overlapping. 

The kinetic theory result (O) has been obtained through the choices Eqns. (0) and (|h and shows no dependence 
on the density. This is a consequence of introducing the factor A~^ in Eqn. ^. We verify in Fig. || that simula- 
tions performed at different number density, with the remaining parameters s, Cl held fixed, show that the thermal 
diffusivity D is indeed independent of n, as predicted in Eqn. (pT|). 



V. CONCLUSIONS 

In this paper we have apphed the model introduced in Ref. pO] to the case of heat conduction in a random sohd. We 
have chosen the equation of state s(ej) for a perfect soUd model and a particular model for Kij, Eqn. (||). It has been 
shown that the model has the correct equilibrium properties and that it reproduces Fourier's law in non-equilibrium 
situations. We have also corroborated, by means of simulations, the kinetic theory predictions for this model which 
are presented elsewhere. Kinetic theory predicts a thermal diffusivity which depends only on the overlapping and not 
on the density. By introducing the factor X~^ in the definition of Kij we have shown that it is possible to interpret n 
not as the density but as the resolution at which the physical problem is being studied. A similar discussion about 
resolution was presented for the case of DPD with no energy conservation in Ref. |jl^ . 
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FIG. 1. The probability density distribution f{e) of the energy of a given particle. Two simulations with q = 50 (left) and 
a — 100 (right) are presented. Other parameters are A'" — 1000, s = 1.5, T^q = 1. Also shown, although not visible because 
they are perfectly superimposed, are the theoretical predictions for the microcanonical and canonical distributions for the given 
values of a. 
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FIG. 2. Average temperature of the particles in each of the 25 bins in which the x axis of the box has been divided. Solid 
line is the theoretical Fourier profile. Open circles correspond to N = 100 and solid circles correspond to A'^ = 1000 particles 
within the box. Other parameters are s = 1.5, Cl ~ 10*. 
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FIG. 3. Average x component of the heat flux as a function of the imposed temperature gradient. A linear dependence is 
obtained with small deviations at large gradients. The slope of the straight line gives the thermal diffusivity. 
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FIG. 4. Thermal diffusivity coefficient for different values of the overlapping. Also shown (solid line) is the kinetic theory 
prediction. 
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FIG. 5. Thermal diffusivity coefficient for different values of the density. Also shown (sohd line) is the kinetic theory 
prediction for this value of the overlapping s = 1.5. 



